Superfluid-insulator transitions of two-species Bosons in an optical lattice 
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We consider a realization of the two-species bosonic Hubbard model with variable interspecies 
interaction and hopping strength. We analyze the superfluid-insulator (SI) transition for the relevant 
parameter regimes and compute the ground state phase diagram for odd filling at commensurate 
densities. We find that in contrast to the even commensurate filling case, the superfluid-insulator 
transition occurs with (a) simultaneous onset of superfluidity of both species or (b) coexistence of 
Mott insulating state of one species and superfluidity of the other or, in the case of unit filling, 
(c) complete depopulation of one species. The superfluid-insulator transition can be first order in 
a large region of the phase diagram. We develop a variational mean-field method which takes into 
account the effect of second order quantum fluctuations on the superfluid-insulator transition and 
corroborate the mean-field phase diagram using a quantum Monte Carlo study. 

PACS numbers: 03.75 Mn, 75.10.Jm, 05.30 Jp 



I. INTRODUCTION 

Experiments with ultracold atoms have achieved re- 
versible tuning of bosonic atoms between superfluid (SF) 
and Mott insulating (MI) states by varying the strength 
of periodic potential produced by standing laser lights. 
The physics of such ultracold atoms in the Mott in- 
sulating state can be described by a bosonic Hubbard 
model, well known in context of other condensed matter 
systems?.. However, ultracold atoms in optical lattices 
offer much better control over microscopic parameters of 
the model. Consequently, it is possible to explore param- 
eter regimes which are not available in other analogous 
condensed matter systems. 

Recently, experiments involving internal states of these 
atoms have been carried outi^. In particular, in Ref. 0, 
the two hyperfine states {\F = 2, uif = — 2) = |1) and 
\F = 1, rriF = —1) = |2)) of 87 Rb atoms have been used 
to create entangled states between atoms in different 
wells of the optical lattice. In these experiments, a 7r/2 
pulse is applied to bosons originally in one of the two 
hyperfine states (say |1)), leaving them in eigenstates of 
a x Q|l) + |2)]/ v / 2), where the a denote Pauli matrices 
corresponding to the two hyperfine states. 

To envisage how such experimental systems are rele- 
vant for realization of a two species Bose-Hubbard model, 
consider an optical lattice created using elliptically polar- 
ized light with polarization angle 0. Since the spin states 
with m s = ±1/2 see potentials V± = Vq sin 2 (fcir± 0), the 



hyperfine states |1) 
given by (see Refs. 



and 1 2) experience potentials Vi(2) 
J|fj|7| for details) 



Vi 



V o sm 2 (kx + 0) 



V 2 = ^-(sm 2 {kx + 0) + 3sm 2 (kx-6)) 



(1) 



Consequently, a change in the polarization angle 8 is 
equivalent to a relative shift of the lattices with respect 



to each other. Since the interaction between the bosons 
is short-ranged, such a shift can be used to control the 
inter-species interaction U'. Note that changing the po- 
larization angle also changes the depth of V2, and there- 
fore the corresponding hopping amplitude t 2 . Hence, sys- 
tems of atoms where state selective optical potentials can 
be implemented may provide ideal test beds for study- 
ing properties of the two species bosonic Hubbard model 
with variable hopping amplitudes and interspecies inter- 
action strength. 

Several theoretical works have discussed realizations 
of novel phases in the two-species system in an optical 
lal lic^^LliLli. Because of the inter-species interaction, 
the Mott phase is divided into regions with different long 
range orders. These phases can be described in terms of 
isospin, a quantum number which describes the occupa- 
tion state of a single site by two components&^iiS. For 
a total occupation 2uq — 1, the states \riQ,nQ — 1) and 
I no — l,n ) correspond to isospin states with S z = 1/2 
and S z — —1/2 respectively. However, at the superfluid 
transition point, which can be approached by decreasing 
the strength of the optical lattice, the isospin description 
breaks down because of strong density fluctuations. The 
isospin quantum number S z , which is given by the differ- 
ence in quantized occupation numbers of the two boson 
species, becomes ill-defined at this point. Nevertheless, 
one can still investigate the effect of such isospin order in 
the Mott state on the superfluid-insulator (SI) transition. 
This is the key issue that we are going to address in this 
work. We note that although there have also been earlier 
studies of the SI transitions from such isospin symmetry 
broken Mott state oV 1 , the phase diagram of the system 
for the entire parameter range has, to the best of our 
knowledge, not been charted out and large parts remain 
to be explored. 

Keeping the above-mentioned experimental and theo- 
retical scenarios in mind, we shall study a two-species 



2 



bosonic Hubbard model described by the Bose-Hubbard 
Hamiltoniar&2ii£ 



u 

~~2 



^ (-t a b\ a bj a + h.cj - fJ,y^ y n ia 

n ia (n ia - 1) + 2A^ nan i2 , (2) 



where a = 1, 2 is the species index, t^) denote hopping 
amplitudes for the two species between nearest neighbor 
sites (ij), the matrix element U denotes on-site intra- 
species Hubbard interaction, U' = XU is the inter-species 
interaction and we have taken the chemical potential /i to 
be the same for both the species. Note that a change in 
optical lattice depth by tuning the laser polarization also 
leads to a relative shift of chemical potential of the two 
species. However this shift is usually small and can al- 
ways be compensated by applying an external magnetic 
field since the two species have different magnetic mo- 
ments. For future convenience, we introduce the ratio 
V — tz/ti an d shall take t 2 <t\ (rj < 1). Our aim is to 
study the different phases of the system as a function of 
/i,A,iiand?7 for odd total filling factor. Also, we shall 
refer to the species index a as the isospin label for the 
bosons with isospin S = 1/2. 

Before proceeding with the analysis, we summarize 
the key results of this work. First, we find that the 
superfluid-insulator transition in systems described by 
Eq. |21 can take place with a) simultaneous onset of su- 
pcrfludity of species 1 and 2 (SF-SF phase) or b) coex- 
istence of Mott insulating phase of species 2 and super- 
fluid phase of species 1 (MI 2 + SFi phase) or c), in the 
case of a unit filling Mott state, depopulation of species 
2 (a-SF phase). Second, for a large region of the phase 
diagram the superfluid-insulator transition occurs with a 
discontinuous jump in the number of each species and is 
therefore first order. Third, there is a second order quan- 
tum phase transition between the a-SF and the SF-SF 
superfluid phases which can be viewed as a no = Mott 
insulator-superfluid transition for the bosons of species 
2. Finally our analysis explicitly demonstrates the ne- 
cessity of including effects of 0(t 2 /U 2 ) quantum fluctu- 
ations (beyond the 0(t/U) mean- field theory) for a cor- 
rect quantitative description of the phase diagram and 
the nature of the phase transitions in the system. 

The organization of the paper is as follows. To put this 
work in perspective, we review the results on the Mott 
phases of the Hamiltonian in Eq. © in Sec. [HI In sec- 
tion lllll we study the SI transition using 0(t/U) mean- 
field theory and also discuss the shortcomings of such a 
theory in the present case. This is followed by Sec. IIVI 
where we implement a canonical transformation method 
which takes into account the effect of quantum fluctua- 
tion to 0(t 2 /U 2 ) on the transition and present a detailed 
phase diagram of the model. This is supplemented by 
quantum Monte Carlo simulations in Section In Sec- 
tion IVII we discuss how the different phases can be de- 
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FIG. 1: Schematic phase diagram of two-species boson model 
in the Mott insulating state for ti = ti = 0. Notice the 
two-fold degeneracy at each site for odd fillings. 



tected experimentally. This is followed by a summary of 
our results in Section lYlII 



II. REVIEW OF MOTT PHASES 

In this section, we review the Mott phases of the two 
species Bose-Hubbard model&Siili. Deep inside the Mott 
phase, for t\ =t 2 = 0, the Hubbard Hamiltonian (Eq. [2J) 
reduces to sum of on-site terms Hi given by 



n ln + — 



2Xniin i2 



(3) 



The phases of Hi are characterized by the ground state 
of the system having an integer number of bosons 
ni, 2 (n/U, A) per site. The phase diagram is shown in 
Fig. ^ Apart from the usual even filling phases where 
n i = n 2i phases with odd filling ni — n 2 = ±1, which has 
no counterpart in single species systems, occur. For the 
rest of this paper, we shall concentrate on phases with 
odd total filling, where each site is doubly degenerate 
(ni — 7i2 = ±1) leading to 2 N degenerate ground states 
for a system with N sites for ti =t 2 = 0. 

At finite hopping strengths, this degeneracy is lifted 
by quantum fluctuations which can be studied by sec- 
ond order perturbation theory. More precisely, one 
can carry out this perturbation theory in the regime 
U,U',\U — U'\ >> ti,t 2 where we are far away from 
both the SU(2) symmetric limit (A — 1) and the vanish- 
ing inter-species interaction limit (A <^ 1). In both these 
limits perturbation theory breaks down. 
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FIG. 2: (Color online) Mean-field phase diagram of the effec- 
tive low energy Hamiltonian in Eq. 0] obtained using pertur- 
bation theory at odd filling for no = 1. The phase diagrams 
for other values of odd no are qualitatively similar in the Mott 
insulating regime. 



To compute the fluctuation correction for a Mott state 
with an odd number no of atoms on each site, we divide 
the system into A and B sublattices (to allow for the 
possibility of an antiferromagnetic phase) and use a trial 
wave- function 



i*)=nn 



IV'i 



ferromagnetic (FM) phase with a ~ 0b = 0, and c) XY 
phases with Q a = ®b =h 0- The nature of the transitions 
between these effective isospin phases can be understood 
by plotting the values of the angles 8a and b across the 
different phase boundaries. These plots are shown in Fig. 

Here the angles 9a and Ob have been plotted for three 
different values of 77. In the XY and FM phases 6 a = 0b 
and the lines are indistinguishable, whereas in the AF 
phase one phase takes on a value tt and the other 0. We 
find that there is always an abrupt jump from the XY 
to the AF phase across the AF phase boundary, suggest- 
ing that AF-XY transition is first order. We do not find 
any canted AF phases. The situation here is analogous 
to the first order melting transition of hard-core bosons 
with next nearest-neighbor interaction at half filling 14 . 
The FM-XY transition, on the other hand, is continuous 
and proceeds via continuous change of a and 0b- 

The phase diagram obtained here agrees qualitatively 
with that of Ref. 0, although there is a quantitative dif- 
ference. In our phase diagram, the tricritical point where 
all the phases meet is at A = 0.25, 77 = instead of A = 
0.5, 77 = 0, as found in Ref. 9. To understand why this dif- 
ference arises, we now map the boson Hamiltonian to an 
effective low-energy spin-model. Defining the isospin op- 
erators Sf = (nu - 7i2i)/2, Sf = (&k&2i + 4A«)/ 2 and 
&i = — ^2i^2i)/2 one obtains an effective XXZ 

model in a magnetic fields 



Hxxz = - [ J -L (Sf Sf + SfS]) + J z SfSj 



-BJ2S* + U(1-\)J2S* 



(5) 



where 

\tPa,b) 



cos ■ 



a, 



|no,no -l)+e^ 



sin — r— I no 



l)«o) 



where |ni, n 2 ) i denotes ri\ and n 2 atoms of species 1 and 2 
at site i. A perturbative calculation yields the 0(t 2 /U 2 ) 
correction to the ground-state energy as a function of the 
angles 0a,b and 4>a,b'- 



Nzt\ 
2U 



(1 + r/ 2 )rio(l + cosOa cosO b ) 



+ (1 - rj 2 )n (cosOA + cos#b) 
+ (1 + f] 2 ) [1 - cos 8 a cos B ] 

+ sm(0 A ) sm(0 B ) cos(0a - <Pb) 



"■o 
2A 

ml 

A 



2- A 



(4) 



where z is the coordination number of the lattice. Mini- 
mizing the Ef with respect to 0a,b and 4>a,b, we obtain 
the phase diagram shown in Fig. [21 for no = 1 as a func- 
tion of the parameters rj and A. This phase diagram 
illustrates presence of three types of phases 9 : a) antifer- 
romagnetic (AF) phase with 0a(b) = 0, 9bia) = ^1 b) 



The exchange couplings Ji , J z and the magnetic field B 
are given by 



J± 



2zt\{\ - r, 2 )n 



U 



At 2 



B 



^ 1 -2A 1 - 



nl-1 



2- A 



(6) 



Note that for no = 1, AF, FM and XY phases meet when 
J z = Jj_ = at A = 0.5,77 = provided one neglects 
the magnetic field term, as done in Ref. 0- However, if 
one retains the magnetic field term, the AF and the FM 
phases will meet when J z + 2B/z = and J± = or 
A = 0.25, 77 = 0. 

The XY phase obtained here is identical to the super- 
fluid counterflow (SCF) phase obtained in Ref. and 
also to the v = 1 bilayer quantum Hall state for small 
layer separation where the layer index plays the role of 
isospin 12,13 . The stiffness energy locking the two order 
parameter phases together in the XY phases can be ob- 
tained from Eq. 0] 



Ps = 



d((j>A - (pi 



t>A=4>E 



NztlrjriQ 
2UX 



sin a sin b 



4 




FIG. 3: Plot of optimum values of 8 a and 9b as a function of 
A for different values of r\. In the XY and FM phases 8 a — Ob 
and the lines for 6a(b) are indistinguishable whereas in the 
AF phase one angle takes on a value tt while the other takes 
on a zero value. We find a discontinuous change in both 6a 
and 6b as the AFM phase is entered, signaling a first-order 
transition. 



Note that the U(l — A) J2i Sgi term in Hxxz is a con- 
stant since = 1/4 for all the states in the low energy 
manifold with S Z i = ±1/2. Hence this term does not 
contribute to the low energy effective Hamiltonian and 
does not play a role in the quantum disordering of the 
XY phase. The disordering of the XY phase due to quan- 
tum fluctuation depends only on the exchange constants 



J z , J|| and B. 



III. 0{t/U) MEAN-FIELD THEORY FOR THE SI 
TRANSITION 

In this section, we shall study the SI transition within 
0(t/U) mean-field theory by constructing a site factor- 
izable variational wavefunction which provides an an- 
alytical albeit qualitative understanding of the transi- 
tion. We shall work in the parameter regime where 
U, U', \U - U'\ > h,t 2 . For the sake of clarity, although 
we shall qualitatively comment on the general case, all 
calculations in this section from here on shall be per- 
formed for two spatial dimensions and no = 1. 

Before carrying out the mean-field analysis, we review 
earlier studies of SI transition for the two species Bosc- 
Hubbard model (Eq. 0). In Ref. HJ 

the SI transition 

has been studied for the case t\ = t 2 but with different 
inter-species and intra-species interaction strengths and 
chemical potentials. This has been done using a standard 



mean-field theory^ corresponding to decoupling the hop- 
ping between sites by introducing order parameter fields 
A Q , i.e. b\ a b ja « b\ a A a + b ja A* a - \ A a \ 2 where the fields 
A a satisfy the self consistency relations A a = (b a ). Their 
analysis led to the prediction of three different phases 
(1) Both species superfluid (SF-SF). (2) Species 1 super- 
fluid and species 2 in a Mott state (SF-MI). (3) Species 
2 superfluid and species 1 in a Mott state. (MI-SF). It 
has been found (erroneously, as we shall see) in Ref. ^3 
that the Mott states are always destabilized by MI-SF or 
SF-MI phases and there is no direct transition from the 
Mott to the SF-SF phase. The transitions are concluded 
to be second order as in the standard single species Bose- 
Hubbard models. The question of the interplay between 
the exchange effects and the SI transition in the region of 
small A was studied by Demler et al& for unit filling fac- 
tor and fixed chemical potential fi/U — hX. Apart from 
the phases mentioned above they also found a superfluid 
phase with species one superfluid and depopulation of 
species 2. 

In our proposed setup, A is not necessarily small and 
the SI transition in this regime has not previously been 
investigated. To analyze the SI transition to 0(t/U) 
within mean field it suffices to consider an on-site trial 
wavefunction 15 

= I[(uo\l,0) i + r\0,l) i +Pi\2,0) i 



-hp 3 |l,l) < +p 3 |0 1 2) i + A 1 |0 1 0) i 



(7) 



where uq is the amplitude of the ferromagnetic Mott state 
1 1,0) in ty v , h and p are amplitudes of removing, adding 
bosons of species 1, 2 to the Mott state and r is the am- 
plitude of isospin-flip. We note here that allowing the 
isospin-ffip process in the trial wavefunction (Eq. is 
absolutely crucial for correctly taking into account the 
manifold of low energy boson states which are degener- 
ate to 0[t/U). The normalization of the wave-function 
yields the constraint Mq + r 2 + p\+ p\+ p 2 + h\ = 1. This 
wavefunction, as we shall see, is appropriate for studying 
the SI transition from the FM and the XY Mott phases. 
We shall comment about the AFM-SI transition later. 

The energy of the variational ground state 
E v (u ,r,p 1 ,hi,p 2 ,h 2 ) = (V V \H\& V ) is given by 



E v = E 



Mott 



£ 

a=l,2 



(pf + p 2 3 )SE+ + p 2 2 SE+ + h\5E^ 



A 2 

a 

ztct 



(8) 



where Emou is the energy of the Mott state, z is the 
coordination number of the lattice, 5E^ denote the en- 
ergies of adding/removing a boson of species a to/from 
the Mott state given by 



6E7 



-u + u, se; 



-u + XU 5E^ = fi, (9) 
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TABLE I: Parameter values of the variational wavefunction 
corresponding to the different phases. 



and the superfluid order parameters A1.2 can be calcu- 
lated from this variational wave-function: 

Ai = zii I = zti (uop-lV2 + rp 2 + u hi) 

A 2 = 2* 2 <*«|6a|*«) = zt 2 (nop* + rhi + rpsVfylO) 

Mathematically, it is possible to show that for all of the 
Mott and superfluid phases (except the AFM phase for 
which we need to use two sublattices), the variational 
energy has a stationary point. The parameter values 
at these points and how they translate into the various 
phases is shown in Table [I] The transition to superfluid- 
ity from the Mott state occurs when it becomes energet- 
ically favorable to have non-zero A a i.e. non-zero ampli- 
tudes of additional particles and holes (p and h) in the 
variational ground state. For our purposes, it is sufficient 
to take all the coefficients real. This amounts to setting 
the phase of the superfluid order parameter to zero and 
does not affect the variational energy of the state. Note 
that the wavefunction (Eq. 01 is general enough to incor- 
porate both the FM (u = 1) and the XY (u = cos(0/2) 
r = sin(6>/2)) phases of the Mott state. However, since 
these two states are degenerate to 0(t/U), our simple 
mean-field treatment cannot distinguish between their 
isospin order. In Sees. IIVI and we shall carry out 
more sophisticated treatments of our model which will 
take into account the effect of 0(t 2 /U 2 ) fluctuations us- 
ing canonical transformation and quantum Monte-Carlo. 
In this section, we shall analyze the 0(t/U) mean-field 
theory and point out certain qualitative features of the 
phase diagram. 



A. General features of the phase diagram 

Although the variational wave function in this section 
excludes second order exchange effects, the qualitative 
features of the SI transition from the FM and the XY 
phases for Uq = 1 can be understood from Eqs. |S] and 
Consider, for example, approaching the SI transition 
from FM/XY Mott phase. For /j < U, XU, since 6E± < 
SE^SE^, at the SI transition point t\ = t\ = SE^/z 
and the energy of the variational wavefunction is mini- 
mized with r = 0, uq r~j 1, pi = P2 = P3 = 0, hi 7^ 0. 
Consequently, from Eq.^J we have A2 = 0, Ai 7^ 0, and 
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FIG. 4: (Color online) 0(t/U) mean-field results for the SI 
transition from the XY phase for A = 0.4 and r\ = 0.8. The 
inset shows 712 at the transition. Superfluidity sets in a) with 
depopulation bosons of species 2 in region A (a-SF phase) 

b) simultaneously for both species in region B (SF-SF) and 

c) with Mott phase for 2 and superfluid for 1 in region C 
(MI2 + SFi). The vertical dotted lines are guides to the eye 
and represent the positions of u c i and /i C 2 (see text). 



also 



Til = (^l&H&H |*«> = 

n 2 = {^ v \b\jb 2 i \®v) = 



-2p\ 



Pi) 



(11) 



Thus the transition to superfluidity occurs with complete 
depopulation of species 2. We refer to this phase as a-SF. 
Alternatively one can view this phase as SF1-MI2 with 
a zero filling factor in the Mott phase. Numerically, we 
find that such a depopulation occurs till a critical value 

of /i = [l cl . 

In the other limit, when [i ~ XU > /i C 2, it is much 
more favorable to destabilize the Mott state by adding a 
particle of species 2 since dE^ <C SE^ , SE^ . As a result 
the transition occurs with uq = 0, r ~ 1 and p 2 7^ 0. 
Consequently, the transition takes place with 



A 2 = 0, n 2 , Ai, m ^ 



(12) 



i.e., we have a transition which is accompanied by a jump 
of population species 2 at the transition. The phase con- 
sists of a Mott insulator of species 2 (since A 2 = 0) and 
superfluid of species 1. We call this state MI 2 + SFi. 

For /i C 2 > /i > fici, SEj and 5E± are comparable 
and the energy of the ground state at the transition is 
minimized for r, uo 7^ and p 2 7^ 0. In this case, provided 
r/ 7^ 0, both Ai and A2 are non-zero at the transition 
implying simultaneous onset of superfluidity of species 1 
and 2 (referred to as SF-SF). The width of this region 
is expected to be large at large 77, since higher t 2 makes 
it energetically more favorable to realize superfluidity of 
species 2. 



6 



The values of fj, c ± and fj, C 2 are shown for representative 
values of 77 = 0.8 and A = 0.4 in Fig. 0] The phase dia- 
gram corroborates the above discussion. From the inset 
of Fig. 21 we find that there are three distinct regions 
where the SI transition takes place with a) depopulation 
of species 2 (region A), b) simultaneous setting of su- 
perfluidity of the two species (region B), and c) Mott 
insulating phase of species 2 and superfluidiy of species 1 
(region C). The situation here is in sharp contrast to the 
even filling case which will always have an intermediate 
state with superfluidity of species 1 and insulating state 
of species 2 for < rj < 1. 

Upon further increase of t\ from the critical value t\ c , 
two scenarios are possible. If superfluidity sets in with 
depopulation, increasing t\ does not change the situation 
further. On the other hand, if the transition occurs to 
either the SF-SF or MI2 + SFi phase, upon increasing 
t\ , the fraction of B atoms in the superfluid decreases as 
shown in Fig. [3] for a set of representative values of 77, A 
and fi. Finally, one crosses a critical value t\ at which it 
becomes energetically favorable for the system to depop- 
ulate. This happens at large enough t\ > t* ~ SE^/z, 
at which the variational energy minima shifts to uq =/= 0, 
r = P2 = 0. Beyond this point, we only find superfluidity 
of species 1. Within 0(t/U) mean-field theory, such a 
transition from SF-SF to a-SF phase is found to be first 
order since ^] a=12 A^/z( a is discontinuous across the 
transition. 

Although we do not show it explicitly here, a similar 
consideration remains valid for the SI transition from the 
AFM phase. This can again be seen by dividing the lat- 
tice into the usual A and B sublattices and constructing 
an appropriate two sublattice variational wavefunction. 
We do not find any translational symmetry broken su- 
perfluid phases. We also note that the above discussions 
have to be modified for no 7^ 1, where the Mott state 
can also be destabilized by adding holes of species 2. For 
example, when no > 1, for fi -C A, 6E^ SE^ ■ Thus 
if 77 7^ 0, we expect the ground state energy to be always 
minimized for Uo,r 5^ at the transition. Consequently, 
there will be no depopulation for any finite 77 in this limit. 
In the rest of this work, we shall restrict all discussion to 
odd fillings with no = 1. 



B. Necessity of going beyond the 0(t/U) mean-field 
theory 

We now discuss the limitations of the mean-field the- 
ory to set the stage for incorporating the fluctuation ef- 
fects. To understand why using the mean-field theory 
is dangerous in the present context, consider plotting 
the mean- field phase diagram at a fixed ji/U = 0.5 and 
ti/U = 0.04 as a function of 77 and A. Such a phase di- 
agram is shown in Fig. [S] Here, we have used the phase 
diagram (Fig. |2J) of the XXZ model [Eq. 0] to deter- 
mine the isospin phases since the 0(t/U) mean-field the- 
ory can not distinguish between them. As can be seen, 
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FIG. 5: (Color online) A plot of the order parameters Ai, A2, 
and 712 for fi/U = 0.57, A = 0.6 and rj = 0.2 as a function of 
the hopping amplitude ti. The system enters the MI2 + SFi 
phase at t — t\ (vertical dotted line between regions A and B) 
from a FM Mott phase (region A). Note that the transition 
occurs with a spontaneous jump of 712 and is hence expected 
to be first order. As we increase ti, both species becomes 
superfluid (region B), until t\ reaches £* (vertical dotted line 
between region B and C) where the system depopulates. The 
depopulation occurs with a jump in Ai and is therefore first 
order. 



the phase diagram (Fig.EJ corroborates the expectations 
based on the qualitative discussion of the Sec. IIII AI For 
small 77 (transition from FM/AFM phases), the system 
favors a-SF phase while for larger 77 (transition from the 
XY-phase) the SF-SF phase dominates. 

However, consider now plotting such a phase diagram 
near \x c \ or /i C 2. Clearly, we expect that incorporating 
exchange effects will make it harder for the isospins to 
flip, since now it costs an energy 0(t 2 /U 2 ). This will, 
in general, shift the positions of ji c i and u C 2 from their 
mean- field values. Therefore, near \x c \ or fi C 2, the phase 
diagrams in the 77 — A plane predicted by the mean-field 
theory will be qualitatively different from the true phase 
diagrams. In this sense, the failure of the 0(t/U) mean- 
field theory in the present case is much more severe com- 
pared to the usual SI transition for single species bosons. 
However, as long as we are away form the critical fi val- 
ues, such a mean-field theory gives qualitatively correct 
results and therefore the scenario described in the previ- 
ous section remains largely valid. 

Another problem of the 0(t/U) mean-field theory is 
that it overestimates the jump of ri\ or n 2 at the transi- 
tion since it does not take into account the energy cost of 
an isospin flip. Consequently, it can erroneously predict 
first-order MI-SF or a-SF to SF-SF transitions where in 
reality such transitions might be second order. Also, as 
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FIG. 7: Bonds in a 2D square lattice. There are two types of 
bonds, horizontal and vertical labeled by and u v respec- 
tively. The horizontal bond shown is denoted by ah- The 
sites on the left and right sides of this bond of (sites i and j 
in this case) are labeled by a hL and a hR respectively. 



0.1 I 1 1 1 1 1 1 1 1 1 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

T| 

FIG. 6: 0(t/U) mean-field phase diagram for the two-species 
model as a function rj and A for fi/U = A/2 and ti/U = 0.04. 
In the absence of second order fluctuation corrections, the SF- 
SF phase borders to the FM phase implying a discontinuous 
change in the population of species 2 for a large region of 
parameter phase. 



we shall see in the next section, the shapes of the tran- 
sition curves and topology of phase boundaries change 
quite a bit upon inclusion of the fluctuation corrections. 

Thus, although the 0(t/U) mean- field theory correctly 
predicts the qualitative nature of the MI-SF transition for 
most parts of the phase diagram it fails drastically either 
when we are close to ji c i or fj, C 2 or when we want to 
estimate the order of the transition. In the next section, 
we remedy this failure by incorporating the 0(t 2 /U 2 ) 
fluctuation corrections. 



IV. CANONICAL TRANSFORMATION 

The effect of fluctuation to second order in t/U can 
be taken into account using a suitable canonical trans- 
formation method. We describe the implementation of 
this method in subsection IIV Al and present the phase 
diagrams in subsection IIV Bl 



A. Implementing the canonical transformation 

We begin by separating the Bose-Hubbard Hamilto- 
nian (Eq.^J into an onsite term Ho = Hi (Eq.^J and 
the hopping terms T. The first step is to write T in terms 
of sum over bonds a between the neighboring lattice site. 
To this end, as shown in Fig. [7| we can decompose the 
hopping into hopping on vertical and horizontal bonds, 
labeled <x V h, between adjacent sites. The hopping term 



can then be written as a sum over bonds 

t = E T -=E T ^+E T -« ( 13 ) 

T<r h = — E *« (C R a^« + kc) (14) 
T °. = -E^( 6 ^<A„ D «+h.c) (15) 

a 

We now seek a unitary transformation that will capture 
the effects of the second order exchange effects. We shall 
only consider the case no = 1, although generalization to 
other values of no is straightforward. To do this, we first 
introduce the projection operators P CT acting on the two 
sites associated with each bond a. P a projects the state 
of the system to the manifold of states having one particle 
on each site of the bond. Such a projection operator can 
be decomposed into two parts depending on whether the 
bosons occupying the sites of the bond are of the same 
or different species: P a = P° + P 1 . For instance for a 
horizontal bond we can write 

Pt„ = P? +Pi 

P° ah = (|i,o)(i,o|) CThi ®(|i,o)<i,o|)^ 

+ (|0,l)(0,l|) CTh ®(|0,l)(0,l|) CTh (16) 

!l L "-R 

Pl h = (|0,l)(0,l|) CThi ®(|l,0)(l,0|)^ 

+ (|l,0)(l,0|) CTh 8(10,1X0,11)^ (17) 

where the state |1,0) denotes |ni = l,7i2 = 0) as before 
and the subscript on the parentheses denotes the bond 
and which of the sites the operators act on (cf. Fig. EJ. 
Using the projection operator (Eas. ITT)il7|) . we now de- 
compose T a (Eq. I13|) into two parts 

T a = PjrT a Pjr + (T a P a + P a T a )=T°+Ti (18) 

where Pj- = 1 — P„. The idea is now to use these results 
to seek a unitary transformation 

H* = e iS He- iS =H + [iS,H} + -[iS,[iS,H}) + ... 

(19) 
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which eliminates the terms to 0(t/U). It turns out 
that a suitable choice of S is 

s = ^E[ AP °+ p - T -] = aZ/E^'^] (20) 

a a 

where we have introduced the notation P£ — AP" + P„ 
for future convenience. We now expand H (Eq. 119(1 in 
powers of S. Since S is first order in t/U , this is equiva- 
lent to an expansion in t/U and we have to 0(t 2 /U 2 ), 



H* = H + Y,T a 

1 



a 

iS, [iS,H ]} + O(t 3 /U 3 ) 



(21) 



We now evaluate the different terms in Eq. (|21jl . 
The algebra is straightforward, but lengthy and we 
present some details in Appendix A. The final result, to 
0(t 2 /U 2 ), is 



H* 



h + J2p«t°p« 

(J 

—2 (P a T a T a+ jP a+ j - T a P a Pv+jTp+j) 



h.c] 

(22) 



where the sum over j extends over bonds which are near- 
est neighbors to a. We note that the third and the 
fourth terms of Eq. [22 represent the effective XXZ model 
[Eq. JSJl] of Section [H] and the two-particle hopping pro- 
cesses respectively, whereas the terms in the last line in- 
volve hopping operators on neighboring bonds and are 
expected to be important in the superfluid phases. 



B. Phase Diagram 

The phase diagram of H* is obtained by dividing the 
lattice into two sublattices A and B and using an on-site 
variational wavefunction \ip v ) = Yl ieA IljeB Wj- 
The division into two sublattices is essential for taking 
into account the AFM phase. We note that this is equiv- 
alent to generalizing the mean-field treatment of Sec. 
IIII to incorporate second order fluctuation corrections. 
Although it is cumbersome to evaluate the expectation 
value {4>v \ H* \tp v ) analytically, it can be calculated nu- 
merically by representing the various operators in the 
Hamiltonian as matrices in an appropriately chosen ba- 
sis. The task of minimizing (ip v \H* \ip v ) is then a nu- 
merical optimization problem. Truncating the Hilbcrt 
space to have at most two particles on each site, we 
perform constrained (to keep the norm to unity) opti- 
mization for each point in the phase diagram. We use 
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FIG. 8: (Color online) Q(t 2 /U 2 



mean-field phase diagram 
for the two-species model as a function 77 and A for fi = XU/2. 
The plot for ti/U = 0.04 should be compared with the phase 
diagram obtained to 0(t/U) in Fig. [(3] (see text). Note the 
gradual evolution of the different phases with increase of t%. 
The existence of the multicritical point is due to the special 
symmetry at /i = \U /2 where adding a hole or a particle of 
species 2 to the FM Mott state cost equal energies. Quantum 
Monte Carlo however reveals that these multicritical points 
can be split (See Sec, fy*!. 



a sequential quadratic programming algorithm from the 
MATLAB (TM) optimization toolbox for this task. Due 
to nontrivial energy landscapes and possible existence of 
first order transitions, several starting points, including 
random starting points, were used as input to the algo- 
rithm. 

First, we show the phase diagram in the 77 — A plane for 
fi/U = 0.5A in Fig. EI which shows the gradual evolution 
of the phases of the system as t\ is increased. A compar- 
ison of Fig. for ti/U = 0.04 to its 0(t/U) mean-field 
counterpart [Fig. 10], immediately shows us the impor- 
tance of incorporating the exchange effects. Whereas the 
0(t/U) mean-field phase diagram shows a large bound- 
ary between the FM and the SF-SF phase indicating a 
first order transition, Fig[H]shows only second order phase 
boundaries and no direct transition between FM and SF- 
SF phases. This clearly points out that incorporating the 
exchange effects can lead to qualitatively different results. 

The transition between the FM and the a-SF phases 
is second order, as expected. The transition between the 
a-SF and SF-SF phases is also found to be second order, 
in contrast to the prediction of the 0(t/U) mean-field 
theory. This is a consequence of incorporating the sec- 
ond order fluctuation corrections. We note that the a-SF 
- SF-SF transition can alternatively be viewed as a Mott- 
insulator (with uq = 0) - superfluid transition of species 
2 in the presence of species 1 in a superfluid state. The 
supersolid (SS) phase obtained for small values of A and 
rj represents a superfluid phase with broken sublattice 
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FIG. 9: (Color online) 0(t 2 /U 2 ) mean-field phase diagram for 
the two-species model as a function r\ and A for (i/U = 0.1A 
and = 0.01. The superfluidity is a-SF for a large param- 
eter regime as predicted by the 0(t/U) mean-field theory. 

symmetry. This is precisely the region where zt 2 /\U be- 
comes large and the perturbation theory breaks down. 
We shall see in the next section using Monte Carlo that 
the SS phase is indeed an artifact and signifies the break- 
down of perturbation theory. 

Similar phase diagrams for fi/U = 0.1A and /i/U = 
0.9A and ti/U = 0.01 are shown in Figs . 151 and 1 1 01 respec- 
tively. These plots confirm that the qualitative expecta- 
tions of the 0(t/U) mean-field theory. For \ijJJ = 0.1A, 
we find a large a-SF region and the transition to a-SF oc- 
curs from both FM and XY phases [Fig. @], the transi- 
tions from XY to a-SF being first order. For fi/U = 0.9A 
(Fig. [TP)), the a-SF phase is replaced by MI 2 + SF 1 . Here 
we find a direct first order transition between the FM 
and MI2 + SFi phases. These first order transitions are 
in perfect agreement with the predictions of the 0(t/U) 
mean- field theory. 



V. QUANTUM MONTE CARLO 

To verify that the inclusion of 0(t 2 /U 2 ) corrections 
using the canonical transformation procedure as carried 
out in section IIVI really gives an improvement over the 
0(t/U) mean-field theory in section ITTT1 we have per- 
formed quantum Monte Carlo (QMC) studies using the 
Stochastic Series Expansion method introduced by Sand- 
wik et ali&ZQ. Here we have used a particular form of 
thes e up dates, directed loop-updates. For details see 
Ref. |2l|. The basis states used were from a truncated 
Hilbert space in which each site hosts at most two atoms 
per site, in the same way as done in the mean-field and 
canonical transformation treatments. We expect such a 
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FIG. 10: (Color online) 0(t 2 /U 2 ) phase diagram for the two- 
species model as a function 77 and A for fj,/U — 0.9A for 
ti/V = 0.02. We see a large area of MI2 + SFi phase again 
in accordance with the mean-field theory prediction. 



truncation to be adequate for reproducing the phase dia- 
gram since we work with no = 1 and are always close 
to the MI-SF transition. We have investigated phase 
diagrams in the range fx/U = 0.1A — 0.8 A and found 
Monte-Carlo results agreeing well with the qualitative 
predictions of both Sees. II I II and IIVI and we now turn 
to a critical comparison between the results obtained by 
different methods. 

First, QMC confirms our suspicion that the appearance 
of the SS phase for small A is indeed an artifact of the 
breakdown of the second order perturbation theory at 
small A. Second, comparisons with QMC show that the 
details of the phase diagram are much better reproduced 
using the canonical transformation method. 

A comparison between the different methods can be 
seen in Fig. II II where the phase diagram has been drawn 
for n/U = 0.25A, h/U = 0.02. The dotted, dashed and 
dash-dotted lines are the phase boundaries as obtained 
using the 0(t/U) mean field theory of section lllll As can 
be seen there is a large discrepancy between the location 
of the phase boundaries as compared to the 0(t 2 /U 2 ) 
mean field theory. The discrete set of points represent 
the phase boundaries obtained by QMC. Comparing the 
three methods we thus see that the inclusion of 0(t 2 /U 2 ) 
effects yields the phase boundaries with great accuracy. 

We also find that the Monte Carlo study predicts the 
a-SF/SF-SF transition to be second order. This vali- 
dates the result obtained using the canonical transfor- 
mation method and shows that the expectation of a first 
order a-SF-SF-SF transition based on 0(t/U) mean-field 
theory (Sec. IHIfl was an artifact of omitting 0(t 2 /U 2 ) 
corrections. 
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FIG. 11: (Color online) Comparative phase diagram for the 
two-species model as a function r\ and A for fj,/U — 0.25A 
for ti/U — 0.02. The discrete points represent the phase 
boundaries as calculated using quantum Monte-Carlo while 
the colored regions are obtained from the 0(t 2 /U 2 ) mean 
field theory. The shape of the markers represent the different 
phase boundaries; a-SF/SF-SF (stars) ; SF-SF/XY (circles); 
a-SF/FM (right triangles); FM/XY (up triangles); a-SF/XY 
(squares). As comparison the phase boundaries obtained us- 
ing the 0(t/U) mean field theory are also shown as lines; a- 
SF/Mott (dotted); a-SF/SF-SF (dash-dotted); SF-SF/Mott 
(dashed). 



An interesting prediction, and possible failure, of the 
0(t 2 /U 2 ) mean- field theory is the existence of continuous 
parameter region having points where four phases meet 
(cf. the diagrams for t\ — 0.04 and t\ — 0.05 in Fig. 0. 
However, using QMC for systems of sizes up to 20 x 20 
sites and inverse temperatures (3 — 1500/Z7 for parameter 
values = 0.5XU, t\ = 0.04, suggests that although the 
phase boundaries come very close they do not meet at 
a single point but a small region showing a first order 
transition between the a-SF and the XY phase seems to 
remain. 



VI. DETECTING THE DIFFERENT PHASES 

The traditional way of examining the existence of su- 
perfluidity in trapped boson systems is to switch off the 
trap, let the cloud of atoms expand freely and image 
the expanding cloud. The momentum distribution of 
the atoms inside the trap can then be inferred by look- 
ing at their position, or equivalently density, distribution 
in the expanded cloud. Since the momentum distribu- 
tion function of the atoms is characterized by the pres- 
ence/absence of coherence peaks in the superfluid/Mott 
insulating states, such a measurement serves as a qual- 



itative probe of the state of the atoms inside the trapi. 
In our proposed setup, however, such a simple expansion 
alone, which can not distinguish between the two species, 
will not be able to distinguish between all the different 
phases. Nevertheless, since the two species have different 
magnetic moments (mj? — —2 and ttif = — 1), it is pos- 
sible to separate them during the expansion using a pair 
of Stern-Gerlach magnets^^. The expanding cloud will 
then be separated into two clouds if both species of atoms 
are present in the system. This, together with the mo- 
mentum distribution measurement, will qualitatively dis- 
tinguish between all the phases obtained in this work. To 
make this statement concrete, let us consider the phase 
diagram shown in Fig. Here the proposed set of mea- 
surements will lead to a) single cloud with no coherence 
peak (FM phase) or b) single cloud with coherence peak 
(a-SF phase), or c) two clouds with no coherence peaks 
(XY phase) or d) two clouds with a coherence peak (SF- 
SF phase). Thus this methods allows, for instance, the 
detection of the mixed phases (MI 2 + SFi, a-SF, and SF- 
SF). It further provides a tool for finding the phase tran- 
sitions between the Mott phases for n = 1. The second 
order transition from the XY-phase to the FM phase, for 
example, will be associated with gradual depletion of the 
atoms in one of the clouds whereas the transition to the 
AFM phase will be characterized by an abrupt change 
from a single to a double cloud. 

The above mentioned detection technique, however, 
does not provide any evidence of the isospin order in 
the Mott states since they can not probe the spatial 
correlations between the atoms at a lattice scale. Such 
correlations can be probed, for example, by tilting the 
optical lattice with a potential gradient 1 . Deep inside 
the Mott phase, such a potential gradient will excite 
the system only if £ = E e ii po ie, where £ is the poten- 
tial energy shift between adjacent lattice due to the field 
gradient and Edi po ie is the dipole formation energ}ii*2&. 
The dipole formation energy will sharply change across 
the phase transition lines between AFM-XY and AFM- 
FM phases and consequently the peak in the excitation 
width, measured in Ref. Q] as a function of the applied 
field gradient, shall show an abrupt shift at the transi- 
tions across these phases. In contrast, there will be a 
gradual shift of the peak position as one moves from the 
XY to the FM phase. Alternatively, isospin order can 
also be measured by probing noise correlation of the ex- 
panding clouds2£i2Si2S. For example, a transition from 
the FM to AFM isospin states will be marked by appear- 
ance of additional peaks in the noise spectrum at half 
the reciprocal lattice vector. The detection of the XY 
phase can also be obtained by the Ramsey spectroscopy 
technique as suggested in Ref. 

Another possible way of detecting the phases is to 
image the expanding cloud by passing a linearly po- 
larized laser beam through it. As shown in Ref. y(J, 
the angle of rotation of the plane of polarization (# ro t) 
of the outgoing laser beam is proportional to the net 
m z along the direction of propagation (x±)of the beam: 
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^rot ~ J dx±m z (x±). ro t can then be easily measured 
by passing the outgoing beam through a crossed polar- 
izer since the intensity of the beam coming out of the 
crossed polarizer is J_ ~ sin 2 (6 I0 t)- We therefore ex- 
pect I- to jump discontinuously across any first order 
transitions such as FM-AFM or XY-ASF phase bound- 
aries and gradually change across second order transi- 
tions such as the FM-a-SF or XY-2SF phase boundaries. 
Of course, such measurements have to be supplemented 
with momentum distribution function measurement to 
distinguish between the superfluid and Mott phases. 

A brief comment on system preparation and validity 
of use of grand canonical ensemble is in order here. One 
should note that although use of a grand canonical en- 
semble in the single species case is a good approximation 
due to the presence of a confining potential (in this case 
the confining potential produces spatial inhomogeneities 
in the filling factor and some of the outer regions will act 
like a particle reservoir for the inner region), this is not 
necessarily the case in the two species system and other 
means may have to be sought. 



VII. CONCLUSIONS 

In conclusion, we have studied the MI-SF transition in 
a consisting of two species of ultracold atoms in an opti- 
cal lattice in a previously unexplored parameter region. 
We have used an 0(t/U) mean-field theory to explain 
the qualitative features of the transition in most regions 
of the phase diagram. This is followed by incorporating 
the 0(t 2 /U 2 ) exchange effects using a canonical transfor- 
mation method and a quantum Monte Carlo calculation. 
All of these methods show that the superfluid-insulator 
transition can occur with either depopulation of species 
2 (a-SF phase) or simultaneous onset of superfluidity of 
both species (SF-SF phase) or Mott insulator of species 2 
coexisting with superfluid of species 1 (MI 2 + SFi phase) 
and can be first order in large regions of the phase dia- 
gram. We have also shown that, whereas some qualita- 
tive features of the SI transition can be obtained from 
0(t/U) mean-field theory, incorporating the 0(t 2 /U 2 ) 
corrections is necessary to deduce the details of the phase 
diagram and order of the transitions between the phases. 
Our quantum Monte Carlo study lends strong support to 
the above-mentioned results and also shows screening of 
bosons of species 2 in the SF-SF phase. We also discussed 
possible experimental tests of some of our predictions. 
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Appendix A 

Here we briefly sketch the derivation of H* (Eq. I22|) 
starting from Eq. To do this we show the detailed 
derivation of two terms Hi = [iS, Ho] and H2 — [iS, T a ]. 
The derivations of all other terms follow in a similar fash- 



A. Hi 

To compute Hi, we use Eq. 1201 to expand S and write 
[iS,H ] = -i-^A^ + P^T.Lffo] 

cr 

(7 

T a P°),H 



(23) 



Now consider the first term in the commutator in Eq. 1231 
Noting that the projection operator P a always projects 
onto the states |1, > or |0, 1 >, we see that we can write 



P=T a H -H P>T a = UP=T a 



(24) 



This is an operator identity guaranteed by the construc- 
tion of the projection operator P a . Other terms in Eq. 
1231 can be written in a similar fashion and we have 

Hi = - J2(P*Tv + P°T a +T a P* + T a P°) 

a 

= - E ( p ° T ° + T ° p °) = - E T l ( 25 ) 

cr a 

where in obtaining the last line we have again used the 
properties of the projection operators P£. Combining 
Eas. mn 1251 and 1201 we get the second term in Eg. 1221 
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B. H 2 

Now we consider the term Hi = [iS, J2 a ^al- For ^ ms ' 
as we shall see, it is useful to define the operator M. a = 
P a T a + T a P a . Then one can write, using Ea.[2fj| 

H2 = -tStEK^I'^vI ( 26 ) 
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Note that unlike H%, here the sum extends over two dif- 
ferent bonds a and a' . Consequently, expansion of Eq. 
EH1 leads to terms which can be classified into two cate- 
gories. The first type of terms involves two hopping oper- 
ators T/j on the same bond while the second involves the 
hopping operators on the different bonds: 

H 2 = H 2a + H 2 b 
o 3 

where the sum over j extend over the bonds which are 
nearest neighbors to a. 

We first consider H 2a . We expand the operators M a 
and P* and use the relation [P a ,P£] = 0. Also, we 
note that all terms of the form P^T a P®T a in such an 
expansion vanish identically. After some straightforward 
algebra, one obtains 

H 2a = \T7 \ p * T l p ° + P a TlP$ - 2T a P*T a ] (29) 



Note that the first two terms of Eq. represent sec- 
ond order virtual hopping processes and thus give the 
t 2 /U terms responsible for the isospin ordering of the 
Mott phases, while the third term represents two particle- 
hopping across a bond with an intermediate virtual state 
of one particle on each side of the bond. 

Next we come to computation of H 2 },. We again ex- 
pand out the operators as before. Here, the crucial iden- 
tity is that any terms of the form P a T a P^ + jT a+ j or 
P a +jT a+ jP£T a vanish as long as a and a + j denotes 
nearest-neighbor bonds. Using this, one gets 

H 2b = -—^2^2 (p*T a T a+j P a+j 

a j 

-T a P*P a+J T a+J + h.c.) (30) 



The other term [iS, [iS, Ho]] in Eq.[2]can be computed 
in a similar fashion. Using all these results, we finally 
obtain Eq. 
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